%% =======================================================================
% Table 2 Top
% =======================================================================
clear all; clc;
currentFile = mfilename( 'fullpath' );
[pathstr,~,~] = fileparts( currentFile );
addpath( fullfile( pathstr, '_individual_data/_prob_densities' ) );
opts = optimset('Display','off');

global cv_1stderror cv_2stderror

cv_1stderror = 1.00;
cv_2stderror = 1.96;

% 1981Q1-1984Q4
data    =  readtable('Individual_PRPGDP.xlsx');

rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>1980); 
kk          = find(data.YEAR<1985); 
nn          = intersect(rr,ll);
nn          = intersect(nn,kk); 
tmp         = [data.PRPGDP7 data.PRPGDP8 data.PRPGDP9 data.PRPGDP10 data.PRPGDP11 data.PRPGDP12];
data_matrix = str2double(tmp);
data_matrix = data_matrix(nn,:);
tmp         = data_matrix;
tmp(isnan(tmp)) = [];
data_matrix = reshape(tmp,length(tmp)/size(data_matrix,2),size(data_matrix,2));


x   = [13:-2.00:3]'; 
mu1  = zeros(size(data_matrix,1),1);
std1 = zeros(size(data_matrix,1),1); 
startingVals = [0 1]';
modelFun =  @(p,x) normpdf(x,p(1),p(2));

for ii=1:size(data_matrix,1)
    y = data_matrix(ii,:)'/100;
    [p, other] = lsqcurvefit(modelFun,startingVals,x,y,[],[],opts);
    mu1(ii)  = p(1);
    std1(ii) = sqrt(p(2));
end

% 1985Q1-1991Q4
data    =  readtable('Individual_PRPGDP.xlsx');

rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>1984); 
kk          = find(data.YEAR<1992); 
nn          = intersect(rr,ll);
nn          = intersect(nn,kk); 
tmp         = [data.PRPGDP7 data.PRPGDP8 data.PRPGDP9 data.PRPGDP10 data.PRPGDP11 data.PRPGDP12];
data_matrix = str2double(tmp);
data_matrix = data_matrix(nn,:);
tmp         = data_matrix;
tmp(isnan(tmp)) = [];
data_matrix = reshape(tmp,length(tmp)/size(data_matrix,2),size(data_matrix,2));

x   = [11:-2.00:1]'; 
mu2 = zeros(size(data_matrix,1),1);
std2= zeros(size(data_matrix,1),1); 
startingVals = [0 1]';
modelFun =  @(p,x) normpdf(x,p(1),p(2));

for ii=1:size(data_matrix,1)
    y = data_matrix(ii,:)'/100;
    [p, other] = lsqcurvefit(modelFun,startingVals,x,y,[],[],opts);
    mu2(ii)  = p(1);
    std2(ii) = sqrt(p(2));
end

% 1992Q1-2013Q4
data    =  readtable('Individual_PRPGDP.xlsx');

rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>1991); 
kk          = find(data.YEAR<2014); 
nn          = intersect(rr,ll);
nn          = intersect(nn,kk); 
tmp         = [data.PRPGDP11 data.PRPGDP12 data.PRPGDP13 data.PRPGDP14 data.PRPGDP15 data.PRPGDP16 data.PRPGDP17...
               data.PRPGDP18 data.PRPGDP19 data.PRPGDP20];
data_matrix = str2double(tmp);
data_matrix = data_matrix(nn,:);
tmp         = data_matrix;
tmp(isnan(tmp)) = [];
data_matrix = reshape(tmp,length(tmp)/size(data_matrix,2),size(data_matrix,2));

x    = [8.45:-1.00:-0.55]'; 
mu3  = zeros(size(data_matrix,1),1);
std3 = zeros(size(data_matrix,1),1); 
startingVals = [0 1]';
modelFun =  @(p,x) normpdf(x,p(1),p(2));

for ii=1:size(data_matrix,1)
    y = data_matrix(ii,:)'/100;
    [p, other] = lsqcurvefit(modelFun,startingVals,x,y,[],[],opts);
    mu3(ii)  = p(1);
    std3(ii) = sqrt(p(2));
end

% 2014Q1-end
data    =  readtable('Individual_PRPGDP.xlsx');

rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>2013); 
nn          = intersect(rr,ll);
tmp         = [data.PRPGDP11 data.PRPGDP12 data.PRPGDP13 data.PRPGDP14 data.PRPGDP15 data.PRPGDP16 data.PRPGDP17...
               data.PRPGDP18 data.PRPGDP19 data.PRPGDP20];
data_matrix = str2double(tmp);
data_matrix = data_matrix(nn,:);
tmp         = data_matrix;
tmp(isnan(tmp)) = [];
data_matrix = reshape(tmp,length(tmp)/size(data_matrix,2),size(data_matrix,2));

x   = [4.20:-0.50:-0.30]'; 
mu4  = zeros(size(data_matrix,1),1);
std4 = zeros(size(data_matrix,1),1); 
startingVals = [0 1]';
modelFun =  @(p,x) normpdf(x,p(1),p(2));

for ii=1:size(data_matrix,1)
    y = data_matrix(ii,:)'/100;
    [p, other] = lsqcurvefit(modelFun,startingVals,x,y,[],[],opts);
    mu4(ii)  = p(1);
    std4(ii) = sqrt(p(2));
end
clc;

mu_all  = [mu1; mu2; mu3; mu4]; 
std_all = [std1; std2; std3; std4]; 

data        = readtable('Individual_PRPGDP.xlsx');
rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>1980); 
nn          = intersect(rr,ll); 
data        = data(nn,:);

A   = [data.YEAR data.QUARTER data.ID str2double(data.PRPGDP11)];
[rows, columns] = find(isnan(A));
unique(rows);
A(rows,:)=[];

cil = mu_all-cv_2stderror.*std_all;
ciu = mu_all+cv_2stderror.*std_all; 
A   = [A(:,1:3) mu_all std_all cil ciu];

data_real = readtable('GDPDEF.xls');
nn        = find(data_real.YEAR>1980);

GDPDEF    = data_real.GDPDEF(nn); 
YEAR      = [1981:1:2018]';

count = zeros(size(A,1),1);
for ii = 1:size(A,1)
    yy      =  A(ii,1)+1;
    tmp     = find(YEAR==yy);
    outcome = GDPDEF(tmp);
    if isempty(outcome) == 1
        break
    end
    if outcome>cil(ii) && outcome<ciu(ii)
        count(ii)=1;
    else
        count(ii)=0;
    end
end

coverage_ratio_2std = sum(count)/length(count);
disp('Share with-in 95% CI bound'); disp(coverage_ratio_2std);

mu_all  = [mu1; mu2; mu3; mu4]; 
std_all = [std1; std2; std3; std4]; 

data        = readtable('Individual_PRPGDP.xlsx');
rr          = find(data.QUARTER==4);
ll          = find(data.YEAR>1980); 
nn          = intersect(rr,ll); 
data        = data(nn,:);

A   = [data.YEAR data.QUARTER data.ID str2double(data.PRPGDP11)];
[rows, columns] = find(isnan(A));
unique(rows);
A(rows,:)=[];

cil = mu_all-cv_1stderror.*std_all;
ciu = mu_all+cv_1stderror.*std_all; 
A   = [A(:,1:3) mu_all std_all cil ciu];

data_real = readtable('GDPDEF.xls');
nn        = find(data_real.YEAR>1980);

GDPDEF    = data_real.GDPDEF(nn); 
YEAR      = [1981:1:2018]';

count = zeros(size(A,1),1);
for ii = 1:size(A,1)
    yy      =  A(ii,1)+1;
    tmp     = find(YEAR==yy);
    outcome = GDPDEF(tmp);
    if isempty(outcome) == 1
        break
    end
    if outcome>cil(ii) && outcome<ciu(ii)
        count(ii)=1;
    else
        count(ii)=0;
    end
end

coverage_ratio_1std = sum(count)/length(count);

clc;
confidence_bounds        = {'SPF-density-implied', 'Giordani-Soderlind-2003'};
Percent95                = coverage_ratio_2std;
Percent66                = coverage_ratio_1std;
table_new                = table(confidence_bounds',[Percent95; 0.72], [Percent66;0.44]); 
writetable(table_new,'_tables/table_2_top.txt','Delimiter',' ')
table_new